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The non-equilibrium dynamics beyond linear response of Weyl semimetals is studied after a sudden 
switching on of a DC electric field. The resulting current is a nonmonotonic function of time, with 
an initial quick increase of polarization current followed by a power-law decay. Particle-hole creation 
a la Schwinger dominates for long times when the conduction current takes over the leading role, 
with the total current increasing again. The conductivity estimated from a dynamical calculation 
within a Drude picture agrees with the one obtained from Kubo’s formula. The full distribution 
function of electron-hole pairs changes from Poissonian for short perturbations to a Gaussian in 
the long perturbation (Landau-Zener) regime. The vacuum persistence probability of high energy 
physics manifests itself in a finite probability of no pair creation and no induced current at all times. 


I. INTRODUCTION 

Condensed matter systems, e.g. graphene, 3D topolog¬ 
ical insulators and Weyl semimetals, provide unique op¬ 
portunity to examine fascinating QED effects, like Klein 
tunneling, Zitterbewegung, chiral anomaly or Schwinger 
pair production, most of which barely accessible to ex¬ 
periment otherwise. In addition to this “fundamental” 
appeal, these phenomena play a crucial role in transport 
properties of these systems. 

Weyl semimetals (WSMs) are 3D materials, which sim¬ 
ilarly to the 2D Dirac electrons in graphene, are char¬ 
acterized by linearly dispersing low energy excitations 
around some points in the Brillouin zonei^— . These Weyl 
points are intersections of nondegenerate bands, and are 
stable against perturbations according to their topolog¬ 
ical nature. The low-energy physics of these materials 
mimic the Weyl fermions well-known from high energy 
physics, giving the name WSM. 

Similarly to clean graphene, when the Fermi energy 
in WSMs is near the Weyl point, there are no charge 
carriers available for transport at zero temperature, since 
the density of states vanishes as ~ close to the Weyl 
point. However, in an applied electric field, particle-hole 
pairs created by the Schwinger mechanism^ contribute to 
transport. 

The non-equilibrium state that evolves after turning 
on an electric electric field can be characterized by the 
statistics of the excitations, and by the induced current. 
As pair creation is described by the Landau-Zener (LZ) 
formula in the strong electric field regime, it is intrinsi¬ 
cally related to the Kibble-Zurek mechanising^— , which 
describes the universal scaling of defect generation in 
driven systems near a critical point. Alas, KZ scaling 
gives only the mean number of excitations, and thus does 
not fully characterize the non-equilibrium state. 

Such a characterization, however, is possible through 


all the higher moments or cumulants, as these contain all 
information about non-local correlations of arbitrary or¬ 
der and entanglement. This is practically equivalent to 
determining the full distribution function of the quan¬ 
tity of interest. Therefore, the full distribution func¬ 
tion of the number of electron-hole pairs is also of in¬ 
terest, yielding additional information about the physics 
of Schwinger pair production. Condensed matter physics 
and cold atomic systems thus provide a unique way to ex¬ 
perimentally detect such quantities^i'i^, beyond the cur¬ 
rent capabilities of high energy physics. These ideas also 
relate to the discipline of full counting statistic o^^'^^ , 
were outstanding experiments measure whole distribu¬ 
tion functionsi^idd, and cumulants up to the 15th order 
e.g. in RefJ^. 

Our results on the time evolution of the current and 
statistics of electron-hole pairs in is summarized in TA¬ 
BLE E The time domain is split into three distinct re¬ 
gions with different behaviour, which we call classical (ul¬ 
trashort), Kubo (short), and Landau-Zener regime (long 
perturbations). 


Time domain 

Classical 

Kubo 

Landau-Zener 



."a « i « \/jeE 

\/ vpeE ^ 

ff pairs (n) 

~ Eh^A 

~ E^t 

~ E^t 

Statistics 

Poissonian 

Poissonian 

Gaussian-like 

Current (j) 

~ EtA^ 

~ E/t 

~ E^t 


TABLE I. The electric field and time dependence of the total 
number of excitations or pairs created (n) and its statistics, 
together with the electric current (j) is shown. A is the mo¬ 
mentum cutoff, E is the electric field. 

The time evolution of the current also allows us to 
conjecture qualitatively the behaviour of the steady 
state current-voltage characteristics. For small volt¬ 
ages, the dynamical calculation combined with Drude 
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theory reproduces the results of Kubo formula calcu¬ 
lations, i.e. the current is proportional to the electric 
field. However, Ohm’s law breaks down for larger volt¬ 
ages and the current-electric field dependence becomes 
non-linear. This critical electric field as well as the non¬ 
linear current-voltage relation are important for possible 
transport experiments in WSMs. 

The paper is structured as follows. First, we introduce 
the model and its solution in section m Then we dis¬ 
cuss the evolution of the current and its implications for 
the steady state conductivity in sections [Hi] and HV] The 
statistics of pair creation is studied in section El and it 
is compared with a complementary measure, the vacuum 
persistence probability, in section ED 


II. ELECTRIC FIELD SWITCH-ON IN A WEYL 
SEMIMETAL 

We consider noninteracting Weyl fermions near a sin¬ 
gle Weyl point. A homogeneous electric field is switched 
on at t = 0, which is described by a time dependent 
vector potential A{t) = (e£'t0(t), 0,0). The time evolu¬ 
tion of a given mode p = (jpx,Py,Pz) is governed by the 
Hamiltonian 


H = Vf[p>-eA{t))-a , (1) 


where cr denotes the vector of Pauli matrices and Vf is the 
Fermi velocity. The spectrum consists of two bands as 

±Vfy^p‘^ +p']_, with p± = -|- pI the perpendicular 

momentum. Initially {t < 0), the system is assumed 
to be in the T = 0 vacuum state, with all modes with 
negative single particle energy filled and positive energy 
modes empty. This effective Weyl theory is valid at low 
energies compared to a high energy cutoff vpA introduced 
for integrals over momentum space whenever necessary. 

At t = 0, the electric field is switched on, and the 
time dependent Schrbdinger equation can be solved ana¬ 
lytically using parabolic cylinder functionsi^^— . The in¬ 
stantaneous eigenenergies form two bands as ±e(p) with 
e{p) = Vf^ijpx - eA{t)Y +p\. 

The current contribution from a given mode p is de¬ 
termined by the mode excitation probability np(t), which 
gives the number of electrons created in the upper band 
due to the electric field and also the holes in the lower 
band, with n-p{t = 0) = 0. The current consists of a con¬ 
duction (intraband) and a polarization (interband) part 
as (jx)p(t) = jp(0 + 


fpit) = -evF 


vpjpx - eEt) 

epit) 


(2np{t) 1) 


jP{t) = 


( 2 ) 

( 3 ) 


The total contribution of a Weyl node is obtained after 
momentum integration. In Eq. the Up independent 
background is discarded, as an empty or fully occupied 


band does not carry curren t^^i^^ . In our non-interacting 
model, the total current, excitation numbers and higher 
cumulants are additive, i.e. given by the sum over the 
Weyl nodes. 

The vanishing gap is a signature of the “criticality” of 
the WSM phase. As such, it exhibits scaling properties, 
which allow us to deduce important properties of the sys¬ 
tem without explicitly solving the Schrbdinger equation. 
The excitation probability of the modes satisfies a scaling 
relation (in units of fi,VF,e. = 1), 

= (4) 

which follows from the time dependent Schrbdinger equa¬ 
tion, and holds for any choice of the dimensionless scaling 
parameter b. The invariants of the scaling transforma¬ 
tion yield the natural dimensionless combinations which 

determine the physics e.g. i = 

etc. The dimensionless time t = -A uniquely classi¬ 
fies the excitation probability as a function of p, where 
tE = ^Jh/vF^E is the time scale related to the electric 
field. Time reversal considerations also give constraint 
on the excitation probabilities^^ 


^p(0 — ^eEt—p(0 t 


( 5 ) 


which means that the excitation probability is symmetric 
with respect to Px = \eEt. Accordingly, in Eq. ®, 
and everywhere where spherical coordinates are used, the 
momentum is measured from (eEt/2,0, 0). That is, p = 
yj{px - eEtj^Y +P±- 

The excitation probability as a function of p is quali¬ 
tatively different in the f ^ 1 and t <C 1 cases (Fig. [1}. 
A perturbative solution valid for t <C 1 isi^ 


Tip — 


{eEhpFf 

iVpP^ 


sin^ 



( 6 ) 


This gives a good approximation for the excitation num¬ 
ber for p ^ eEt. At short times high energy states may 
become excited, which is reflected in the power law de¬ 
cay of excitations as a function of momentum (^ p~^ for 
p <C h/vpi). 

If the perturbation is long, the probability of exciting 
a given mode is well approximated by the LZ solutioit^^. 


Up 


<d{px)<d{eEt - Px) exp 


( FVfp]_ \ 

\ heE ) 


( 7 ) 


This describes a “dynamical steady state”, which is char¬ 
acterized by a longitudinally growing cylinder of excited 

states of length eEt and radius ^ In contrast to 

the short time limit, the excitation probability decays ex¬ 
ponentially for large momentum. This exponential decay 
can be explained as a tunneling effect within the WKB 
approach^^. 

Along with the analytical calculations, for comparison, 
we determine numerically Up and dfUp by applying an 
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explicit Runge-Kutta method^i to solve the time depen¬ 
dent Schrddinger equation. In Fig. [1] we compare the ap¬ 
proximations used for Up with the numerically obtained 
values. 


numerical (t = 0.01) 



-5 0 5 


Px/eEt 



FIG. 1. Illustration of the excitation probabilities for short 
(left) and long perturbations (right). The excitation map has 
a “dipolar” character for short times, and the approximate 
formula m is nearly indistinguishable from the numerical so¬ 
lution for p eEt. The excitation map is cylindrical for 
long times. An (asymptotically irrelevant) increased number 
of excitations at pa, = 0 and p^ = eEt is not captured in the 
approximation 0. 


III. EVOLUTION OF THE CURRENT 


We are now in a position to discuss the time evolution 
of the current. The high energy cutoff, up A defines an 
ultrashort timescale tA = which satisfies ^ 
for both condensed matter— and cold atomioi^ realiza¬ 
tions of WSMs, similarly to the 2D case^. The scaling 
property (0 implies a scaling for the integrated current 
as 




( 8 ) 


The particular choice of 6 = Ie allows us to reveal 
the electric field and time dependence of the physi¬ 
cal quantities. The current is expressed as jE^it) = 
^ 3 / 2 jc/p The scaling functions Jify (a;) are de¬ 

termined from Eqs. (1^ after evaluating the momentum 
integrals with the particular form of np(t), 

r-(^)"ln^ t«tA 

Jb.aW (9) 

tE<^t 


J£;,a(0 


E3/2 


t i Ej 

ti 


(1 -|- non-univ.)^ 
const 


t <^t A 

tA < i (10) 

tE 


The term ”non-univ.” in the second line of Eq. (TTO)) de¬ 
notes the non-universal contribution from the high en¬ 
ergy regularization, which dies out with increasing time, 
as discussed further in Eq. (na. 

For t <C is the current is dominated by the polariza¬ 
tion part. Because of the large weight of high energy 
states available to excite at ultrashort times t < Ia, the 
current is determined by the cutoff. The ultrashort time 
contribution of a Weyl point to the current is linear in 
time. 


This behavior has also been observed for 2D Dirac 
fermionsiS, and can be explained using a classical picture 

d^€ 

of particles with effective mass m~j = Qp.^. accelerat¬ 
ing in an external field satisfying Newton’s equation. In 
2D, the current saturates at f ~ ^a, and remains constant 
until t tE- In 3D the behavior is sharply different as 
the current starts to decay as after reaching a max¬ 
imal value at t tA- The precise form of the decay 
depends on the microscopic details (i.e. on the cutoff), 
but the exponent is a universal characteristic of Weyl 
physics. Imposing a sharp cutoff results in an oscillat¬ 
ing current j ~ U^(l -|- cos(t/tA)), also obtained within 
linear response^^. A smooth (exponential or Gaussian) 
cutoff of the form exp(—v^p/A) or exp(—p^/A^) does 
not generate the oscillating part, and gives 

= (12) 

DTT^ nypt 


where F{x) ~ for a; <C 1 and F{x) = 1/2 for x » 1. 
The qualitative difference between the 2D and 3D cases 
is a consequence of their respective phase space sizes. 
The polarization current is a sum of contributions with 
oscillating signs j ~ f dp sin(2pt)p‘^“^, which, by simple 
scaling, results in a time independent contribution in 2D, 
but decays as in 3D. 

The conduction part overtakes the polarization term 
at t ~ tE, beyond which the current increases linearly 
with time and nonlinearly with electric field as 


J(t) 


1 e^F^ 


(13) 


This is simply the number of charge carriers per unit 
volume in the steady-state cylinder multiplied by evF- 
It is beyond linear response, as it depends quadratically 
on the external field^S. Our analytical predictions for 
the current are illustrated on Fig. [2l together with the 
numerical results. 

Bloch oscillations appear on a timescale teioch ~ 
where a is the lattice constant, and our description holds 
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jm 



ln{t) 


FIG. 2. Time evolution of the total current after switching 
on an electric field. The analytical curve is the sum of po¬ 
larization current (HI. dominant for t <^tE, and conduction 
current US, dominant for t. The evolution of the number of 
pairs, Ki, is shown in Fig. [3] 


for t fBloch- The timescale related to the cutoff is 
non-universal, and both Ie and teioch depend on the ap¬ 
plied field. These three scales are in fact not indepen¬ 
dent, which can be seen in the following way. The mo¬ 
mentum cutoff is proportional to the largest momentum 
in the system A = —, which relates the timescales as 
^A^Bloch = c4. where c > 1 is a system dependent con¬ 
stant describing the ratio of the linear size of the Bril- 
louin zone and the validity range of Weyl physics. This 
also implies that in the experimentally relevant t^^tE 
case, Ie «C ^Bioch is also satisfied, and all three regions 
appear before Bloch oscillations set in. It is interesting 
to note that the maximal current is jmax e^VEEA/h?, 
which the system reaches upon leaving the classical re¬ 
gion during the time evolution. Even in the nonlinear 
region in Eq. m, the maximal current falls to the same 
order of magnitude, which is in sharp contrast to 2D 
Dirac semimetals, where the non-linear current strongly 
exceeds the current from the classical region. 

As the external field induces current, it also injects 
energy into the system. The conduction and the polar¬ 
ization current decompose the total pumped energy into 
reversible (“work”) and irreversible (“heat”), as follows. 
An infinitesimal change in the energy can be written as 
dE = + Cidrii), where i = (p, ±) runs over 

all momenta of the two bands. The first term corre¬ 
sponds to the reversible work done on the system, dW = 
'^pdt€p{2np — l)dt = VEj^{t)dt, while the second cor¬ 
responds to the heat exchange, dQ = 2 epdtripdt = 
VEj^{t)dt, where we have expressed everything by the 
properties of the lower band. Correspondingly the work 
done on the system and the heat are 


W = VE 

[ dsf{s) 

10 

(14) 

Q = VE 

/ dsjP(s). 

/o 

(15) 


It is easy to check that the sum of the heat and work 
yields the total energy of the time evolved state AE = 


W -\- Q = ^p 2epnp, i.e. simply the sum of the energy 
absorbed by the excited modes. 


IV. STEADY STATE PICTURE FROM DRUDE 
THEORY 


The Drude picture provides a heuristic way to relate 
our results to optical conductivity studies of a WSM in 
the presence of impurities. In general, this is expected 
to worb^i for (contributions to) quantities independent of 
the relaxation time, as e.g. the high frequency conductiv¬ 
ity or the universal minimal conductivity of graphene^^. 
In this spirit, the dynamics described above holds until a 
characteristic time determined by an effective scattering 
rate 1/t, and the zero frequency limit of the AC conduc¬ 
tivity can be estimated by substituting time as f —>■ r. 
This results in the counterintuitive conclusion that in the 
< T < Ie region, the conductivity is proportional to 
the scattering rate, a{uj —>• 0 ) ~ ’ '''^hich agrees 

with the results of Ref.— based on Kubo formula calcula¬ 
tions. Although this simple Drude picture works well for 
graphenei^, it fails to describe the transport properties 
of WSMs because in 3D, the density of states at the Weyl 
point vanishes even in the presence of small amounts of 
disorder—, and concomitantly the quasiparticle lifetime 
diverges2i^. The Drude picture can be rescued if we apply 
it to jp, substituting the time variable with the momen¬ 
tum dependent lifetime, and then evaluating the integral. 
The scattering rate in the Boltzmann or Born approxi¬ 
mation is 1/Tp = 2TT'yg{ep)^, where g{e) = is 

the density of states, and 7 characterizes the scattering 
strength. In the large scattering limit 7 » , integrat¬ 

ing Eq. ([3]) with np(Tp) from (| 6 ]) reproduces the results 
of Refsi^i^, that is 


cr 


2„,2 


e V 


Ivy 


(16) 


with a different prefactor and an additional logarithmic 
correction ~ ln(^)^). The above treatment is valid 

for small electric fields eE <C , when the dominant 

ri V p 

contribution to the current comes from the momenta sat¬ 
isfying Tp < tfi. 

If the scattering strength is small, such that there is 
enough time for the modes to go through the LZ tran¬ 
sition, then the steady state occupation profile will be 
qualitatively similar to the LZ solution. As the quasi¬ 
particle lifetime is finite everywhere except in the close 
vicinity of the Weyl point, the cylinder of densely excited 
states will not extend to infinity, but will be characterized 
by a finite length eETes{E). The precise form of Teg de¬ 
pends on the detailed nature of the scattering process. If 
there is a constant scattering rate I/r, then Tgg = r, but 
generally it will depend on the electric field. The Drude 
picture estimates the stationary current in the non-linear 
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regime as 


Jstac — 


1 

47r3 


es{E), 


(17) 


and Ohm’s law breaks down. The explicit E dependence, 
however, depends strongly on the precise form of Tf.Q[E). 
In case the relaxation time becomes independent of the 
electric field in the non-linear region, a crossover from 
the j ^ E linear region to a j ^ E"^ non-linear region is 
expected. 


V. STATISTICS OF PAIR CREATION 


The expectation value and time evolution of the cur¬ 
rent is largely influenced by the number of pairs created, 
as follows from Eqs. @, ®. This we now investigate in 
more detail. Although the expectation value of a quan¬ 
tity reveals much about underlying physical processes, 
fluctuations contain essential information as well and are 
important to provide a comprehensive description of the 
system^E. Therefore, beyond simple expectation values, 
we study the fluctuations of the pairs created by their full 
distribution function. This provides a complementary 
measure to characterize the out-of-equilibrium state. As 
opposed to calculating the probability distribution func¬ 
tion of pairs created directly, it is more convenient to 
work with the cumulant generating function (CGF) in 
unit volume, which is the logarithm of the characteristic 
function = y In ^exp(i<piV)y Here, N denotes the 

excitation number operator, and the expectation value is 
taken with the time evolved initial state. The CGF is 
expressed as sum over momentum space, 

7 X! + (exp(i(/j) - l)np] (18) 

p 


The probability density function is the inverse Fourier 
transform of the characteristic function, that is, p{n) = 
^ J dip e:Kp{V4>{ip) — imp). For short perturbation, i.e. 
t tE the excitations add up from an extended region in 
momentum space with small excitation probability. The 
contribution from p < 2eEt, where rip ~ 1, is negligible 
because of the small volume of the domain ^ and 
a Taylor expansion of the logarithm in Fq. (HI gives a 
good approximation, 4>{<p) = {exp{i(p) — 1)'^X)''^P' That 
is, the total number of excitations per unit volume is 
Poissonian as p{n) = A" exp(—A)/n! with 


A = 


1 {eE)'^t 

127 r ^ E^vf 




(19) 


where S^iy) = Jq dx sin^x/x^ = ?/ for ?/ <C 1, while 
it saturates to 7r/2 for i/ ^ 1. All cumulants of the 
Poisson distribution are equal to the single parameter A. 
The first cumulant is the expectation value, that is, for 
t «C tA the excitations are created quadratically in time, 
while for Ia ‘C t -C t_E, the creation rate is constant. This 


behavior is clearly seen in Fig. |3l where we compare the 
numerically determined cumulants with the approximate 
solutions. 





ln(t) 


FIG. 3. Time evolution of the cumulants (ki_ 4 ) of electron- 
hole pairs per unit volume (log-log plot, numerical results). 
The cumulants coincide for t <^tE, which is a clear signature 
of a Poissonian distribution. The grey dashed lines show the 
t ^ tA and t ^ Ia asymptotics of the analytical formula 
Eq. (I19II . For long times, the cumulants branch and follow the 
approximate formulae derived from Eq. (I20II (colored dashed 
lines) within a time independent constant coming from the 
difference between the exact rip and the LZ approximation 

(Fig-in. 


For Ie the excited states are confined to a cylinder 
in momentum space, and substituting Eq. 0 into dUD 
yields 


(f){p) = —aLi2(l — exp(iv3)) 
_ 1 {eE)H 
E^Vf 


( 20 ) 

( 21 ) 


where Li 2 (a:) = dilogarithm 

function^, in agreement with Ref^i. As time evolves 
the higher cumulants start to deviate from the first one, 
and the distribution is no longer Poissonian (see Fig.®. 
The cumulants are determined from the series expansion 
of the CGF, the first few being ki = a, K 2 = a/2, 
^3 = a/6, Ki = 0. Except for the variance all even 
cumulants vanish. There is a time independent contri¬ 
bution from the Px ~ d and px ~ eEt regions in rip 
(see Fig. ®, which is not captured in Eq. 0, which 
gets overwhelmed by the time dependent terms. Apart 
from this, the cumulants listed above approximate very 
well the numerical results (Fig. [3]). The peak of the dis¬ 
tribution function is well captured in the central limit 
theorem (CLT) approximation, which states that the ex¬ 
citation number is Gaussian with mean a and variance 
cr^ = a/2: p{n) = --jl=exp{—{n—a)'^/a}. This ap¬ 
proximation neglects the cumulants higher than the sec¬ 
ond. The asymptotic decay of the distribution can be 
determined from the Gartner-Ellis theorem^, which in 
this case is essentially a saddle point approximation of 
the inverse Fourier transform of the characteristic func¬ 
tion. The probability of exciting a large number of pairs 
decays slower than estimated from the CLT, but still in a 
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Gaussian manner p{n) ^ exp{—n^/2a} (note the factor 
2 difference in the denominator of the exponential with 
respect to the Gaussian distribution). 


VI. PROBABILITY OF NO CURRENT AND 
THE VACUUM PERSISTENCE PROBABILITY 

In spite of the applied electric field, there is a finite 
probability of no pair creation and no induced current, 
also known as the vacuum persistence probability. Ana¬ 
lyzing the decay of this probability provides an alterna¬ 
tive way to determine the pair-creation rate, which was 
used e.g. by Schwinger in his seminal paper—. The vac¬ 
uum persistence probability is Po{t) = | (6| U{t,0) |0) p, 
where U(t, 0) is the time evolution operator in the ex¬ 
ternal field, |0) and |6) are the (Schrodinger) vacua at 
time 0 and t respectively. With the knowledge of rip it is 
expressed as 


Pq = exp 1^— ^ ln(l — np)j = exp {—Vwt ), (22) 

where 

T I A fort<t£, 

(23) 


— T ^ an 


for tE 


is the rate of vacuum decay per unit volume and time, in¬ 
creasing as E'^ and being independent of time for t ^ tA 
and increasing linearly with time for t t a- Alterna¬ 
tively, the pair-creation rate can also be defined as the 
total number of pairs created divided by the time it took, 
i.e. as Ki/t. Nevertheless, these two definitions agree in 
the short time limit and only differ by a constant prefac¬ 
tor in the long perturbation limit (Fig. 2]). The vacuum 
persistence probability characterizes the time evolution 
similarly to the Loschmidt echo^: it measures the over¬ 
lap of the non-equilibrium time evolved wave function 
U{t,0) |0) with a reference wave function, which in this 
case is the adiabatically evolved state. 

So far we have assumed the initial state to be the 
ground state without any excitations, which describes the 
zero temperature response of WSMs. An arbitrary initial 
distribution function can be handled similarly, as long the 
modes with different momenta are independent, which is 
the case e.g. at finite temperature. Let /(p) be the 
probability distribution function of having an excitation 
with momentum p in the initial state. The post-quench 
occupation number is expressed as a weighted sum of 
the excitation probability of an unexcited and an excited 
mode asi^ 


= [1 - /(P - eEt)] rip -h /(p - eEf) [1 - iip] , 

where /(p) = l/(exp(/3ep) + 1), fi = l/fc^T. The initial 
number of excited states due to thermal fluctuations is 
riT ^ l/{f3hvFY, which is small near T = 0, and does 


Rate(jA) 



ln(t) 


FIG. 4. Particle creation rate as a function of time estimated 
from the vacuum persistence probability and from the total 
number of excitations (log-log plot). The dashed lines show 
the results of Eqs. (I19I21II and Eq. (I23II . 


not modify qualitatively the results. This argument ap¬ 
plies for systems with thermal initial density matrices, 
which are detached from the environment during time 
evolution. This assumption needs a thermalization time 
much longer than the observation time, which is usually 
not satisfied in condensed matter, but could be achieved 
with cold atoms. Similarly a small deviation in the Fermi 
energy from the Weyl point gives only a subleading cor¬ 
rection. 


VII. CONCLUSION 

We have studied the nonlinear response of WSMs af¬ 
ter switching on an external electric field before Bloch 
oscillations set in. The ultrashort time dynamics is 
non-universal and the current depends on the details of 
the band structure at high energies. The current and 
the number of created electron-hole pairs grow linearly 
and quadratically with time, respectively. The univer¬ 
sal properties of Weyl nodes are manifested in the in¬ 
termediate and long time responses. In particular, at 
intermediate times, the current decays as 1/t due to the 
interplay of the number of created pairs and the avail¬ 
able phase space. Particles are created at a constant 
rate, generating a Poissonian distribution for the statis¬ 
tics of the electron-hole pairs. At long times, the particle 
creation rate takes on a constant value again, but the 
current starts to increase again in time because of the 
increasingly large number of pairs moving in the same di¬ 
rection. The distribution function of excitations crosses 
over from a Poissonian profile to a Gaussian distribution, 
which follows from the central limit theorem, applicable 
in the long time limit due to the large number of pairs 
created. The real time evolution of the current is trans¬ 
lated to the conductivity of disordered WSMs within a 
generalized Drude picture, reproducing the results of pre¬ 
vious calculations with different methods. This is a re¬ 
markable example of a problem from high energy physics 
which naturally corresponds to one in condensed mat- 
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ter physics with a separate set of observables, and which 
allows an exquisitely detailed analysis, thus holding the 
promise of a detailed experimental study in a tabletop 
experiment. 
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